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ABSTRACT 

The  work  reported  here  Is  part  of  a  program  [3]  for  the 
development  of  mathematical  methods  for  calculations  in  fluid 
dynamics.   In  part  A,  the  work  to  date  on  the  Mach  reflection 
calculation  is  summarized.   In  the  flow  there  are  two  curved 
shocks^  a  slip  surface,  a  plane  shock,  and  rigid  walls.   The 
desired  solution  of  the  fluid-dynamics  equations  is  stationary 
in  similarity  variables;  approximate  initial  data  are  therefore 
assumed,  with  the  expectation  that  the  flow  will  settle  down  to 
the  desired  one  asymptotically.   The  Eulerian  equations  for  the 
smooth  parts  of  the  flow  are  coupled  to  the  jump  and  boundary 
conditions  on  the  various  surfaces.   No  pseudo-viscosity  is  used; 
instead,  two-dimensional  fitting  procedures  have  been  devised  to 
apply  the  jump  and  boundary  conditions.   All  equations  are  accurate 
to  the  second  order  in  the  space  intervals.   Numerous  changes  have 
been  made  during  the  year  since  the  original  formulation  was  coded 
and  debugged:  l)  a  least-squares  subroutine  used  for  quadratic 
fitting  of  functions  of  x   and  y  had  to  be  modified  to  prevent 
loss  of  significant  digits;  2)  it  was  found  necessary  (or  at  least 
highly  desirable)  to  improve  the  initial  data  so  that  all  jump 
conditions  were  satisfied  at   t  =  0  (this  required,  among  other 
things,  constructing  subroutines  that  compute  three-shock 
configurations);  3)  it  was  found  necessary  to  allow  the  number 
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of  points  on  the  slip  surface  to  Increase  during  the  run; 
k)    to  stabilize  the  numerical  treatment  of  the  shocks.  It  was 
necessary  to  add  terms  to  the  Ranklne-Hugonlot  equations  corre- 
sponding to  a  surface  tension  (artificial,  of  course),  which, 
however,  contains  the  square  of  the  separation  of  points  and 
therefore  does  not  Interfere  with  second  order  accuracy; 
'5)  similar  terms  had  to  be  added  for  the  slip  surface  -  In 
this  case,  the  added  terms  do  not  vanish  as  the  separation  goes 
to  zero,  because  they  are  needed  to  stabilize  against  Helmholtz 
Instability,  which  Is  a  physical,  not  a  mathematical  phenomenon; 
6)  many  changes  had  to  be  made  In  the  recipes  for  choosing  point 
clusters  for  the  fitting.   With  these  changes,  all  parts  of  the 
code  now  appear  to  operate  satisfactorily  (although  there  Is  no 
estimate  of  accuracy,  yet)  except  for  the  treatment  of  the  slip 
surface,  which  showed  Instabilities  after  about  50-70  cycles. 
The  main  problem  was  temporarily  discontinued.  In  order  to 
concentrate  on  the  slip  surface.   Part  B  of  this  report  describes 
a  simplified  problem  of  a  nearly  plane  slip  surface  between  parallel 
walls  and  with  a  periodicity  condition  along  the  surface.   This 
problem  uses  all  the  machinery  of  the  main  code  except  for  the 
parts  having  to  do  with  the  shocks,  and  It  has  additional  subroutines 
for  computing  an  analytic  solution  based  on  the  linearized  theory  of 
Helmholtz  Instability  modified  by  compressibility,  surface  tension, 
and  the  presence  of  the  walls.   These  subroutines  are  used  to 
provide  Initial  data  and  automatic  subsequent  comparison  (printed 
out)  of  the  numerical  and  analytic  solutions.   Illustrative 
analytic  results  are  presented.   Complex  arithmetic  Is  used  In 
the  analytic  work.   It  Is  hoped  that  the  experience  to  be  gained 
with  the  simplified  problem  will  show  how  the  treatment  of  the 
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slip  surface  must  be  modified  in  the  Mach  reflection  problem. 
Part   C   of  this  report  indicates  a  number  of  ways  in  which  it 
is  planned  to  improve  the  overall  formulation  at  the  time  of 
conversion  to  or  receding  for  the  7090. 
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PART  A:  Summary  of  Work  to  Date  on  the  MACH  code . 

1 .  Introduction. 

We  consider  a  plane  shock  moving  down  a  shock  tube  as  at 
the  left  of  Figure  1.   We  neglect  viscosity  and  other  dlsslpatlve 
mechanisms,  and  we  assume  a  perfect-gas  equation  of  state.   After 
the  shock  has  struck  a  wedge  or  ramp,  there  are,  under  some  circum- 
stances, a  curved  reflected  shock,  a  Mach  shock,  which  may  be  re- 
garded as  the  partial  coalescence  of  the  primary  and  reflected 
shocks,  and  a  slip  surface  (sometimes  called  a  "vortex  sheet"  ), 
as  at  the  right  of  Figure  1.   For  discussion  of  the  experimental 
and  theoretical  aspects  of  these  phenomena,  see  [l]  and  [2]. 

If  X,  Y  are  cartesian  coordinates  in  the  plane  of  Figure  1, 
with  the  origin  at  the  tip  of  the  wedge,  and   t   is  time  measured 
from  the  Instant  at  which  the  primary  shock  touches  this  tip,  the 
flow  pattern  is  stationary  in  similarity  variables   X/t,   Y/t, 
(until  the  reflected  shock  strikes  the  upper  wall  of  the  shock 
tube).   It  is  therefore  assumed  that  the  initial -value  problem 
of  this  flow,  when  expressed  in  these  variables,  has  the  property 
that  with  almost  any  initial  data  the  solution  approaches  the 
desired  one  as   t  — >  0°.   Physical  considerations  suggest  that 
if  the  initial  flow  pattern  is  not  too  far  from  the  correct  one 
this  approach  will  be  fairly  rapid. 

This  problem  was  selected  for  study  in  connection  with  a 
program  of  research  [3]  whose  aim  is  the  development  of  computing 
methods  for  fluid  dynamics  In  two  space  variables  and  time.   The 
general  ideas  of  the  program  are: 

(l)  When  there  are  more  space  variables  than  one,  the 
Eulerlan  formulation  should  be  used,  even  though  it 
Is  less  convenient  for  the  treatment  of  Interfaces  and 
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contact  discontinuities,  because  the  distortion  of  a 
Lagrangean  point  net  makes  the  difference  equations 
Inaccurate . 

(2)  By  careful  centering,  etc.,  the  difference  equations 
should  be  made  accurate  to  second  order  of  the  space 
and  time  Intervals.  Lax  and  Wendroff  have  shown  [^] 
that  stable  difference  equations  of  this  kind  can  be 
constructed  for  the  Eulerlan  as  well  as  the  Lagrangean 
formulation. 

(3)  The  treatment  of  jump  conditions  and  boundary  conditions 
must  also  be  accurate  to  second  order.   Each  moving 
boundary  must  be  followed,  by  somehow  specifying  its 
shape  and  position  as  functions  of  time,  and  the  jump 
conditions  must  be  applied  by  accurate  fitting  procedures. 

One  of  the  main  tasks  of  the  program  is  the  development  of 
the  fitting  procedures.   The  strategy  is  in  a  sense  opposite  to 
that  of  von  Neumann  and  Richtmyer  [5,6],  who  Introduced  a  pseudo- 
viscosity  so  as  to  smear  out  the  shocks  slightly  but  treat  them 
easily  and  automatically  at  the  expense  of  some  loss  of  accuracy. 
For  problems  in  one  space  variable,  one  can  afford  to  use  a  suf- 
ficiently fine  point  net  that  the  loss  of  accuracy  is  acceptable, 
but  in  two  or  more  space  variables,  the  loss  Is  much  more  serious. 
On  the  other  hand,  the  flexibility  and  capacity  of  today's  computers 
make  it  easy  to  handle  the  rather  complicated  and  sophisticated 
formulas  needed  for  the  fitting  procedures. 
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2.  The  partial  differential  equations . 

In  the  present  problem,  the  equations  for  the  smooth  part 
of  the  flow  are  written  In  terms  of  the  dimenslonless  similarity 
variables 


u  =  u(x,y,t)  =  g-  ,  V  =  v(x,y,t)  =  g-  , 


(1) 


where  X,  Y,   and  U,  V  are  cartesian  coordinates  and  velocity 
components,  respectively.   S   Is  a  constant  having  the  dimensions 


of  velocity  -  for  example,  the  speed  of  the  primary  shock.   We 

dx    bj 


use  vector  notation  x  =  (x,y),  v  =  (u,v),  V  =  ( —  ,  — ),   and 


we  denote  by  D   the  operator 


D^  =  A  +  (u-x)  -^  +  (v-y)  i-  =  A  +  (v-x).V  (2) 

St         Sx         dy    St    ^  "^ 


which  gives  the  time  derivative  along  fluid  particle  paths.   The 
equations  are 


D^  V  +  -\   Vp  =  0,  (3) 

D  p  +  p  V.v  =  0,  (4) 

°o^-  ^°o  P  =  °^  ^5) 

p 

where   p,  p,   and   C.   are  the  pressure,  density,  and  specific 

energy  they  satisfy  an  equation  of  state,  which  we  take  to 

be  simply 
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(7-l)p 


(7  =  constant) .  .  (6) 


5.  Difference  equations  for  regular  net  points. 

In  the  region  between  the  reflected  shock  and  the  slip 
surface,  and  In  the  region  between  the  slip  surface  and  the  Mach 
shock,  a  rectangular  net  of  points  Is  used  for  the  computation. 
For  a  regular  net  point  point  (one  whose  four  nearest  neighbors 
are  In  the  same  region),  centered  difference  quotients  are  used 
for  spatial  derivatives. 

The  difference  expressions  that  were  chosen  to  approximate 
the  left  members  of  (3),  (^),  and  (5)   are  accurate  to  second 
order  in  the  space  differences  Ax  and  Ay  but  only  to  the 
first  order  in  the  time  difference  At.   This  was  considered 
adequate,  since  we  are  Interested  only  in  the  asymptotic  stationary 
solution.   (On  the  other  hand.  Lax  and  Wendroff  have  shown  [^] 
that  the  equations  can  be  made  accurate  also  to  second  order  in 
At;   I  now  feel  that  the  choice  made  in  the  MACH  code  was  somewhat 
unwise,  since  the  Lax -Wendroff  equations  are  not  much  more  complicated 
than  the  ones  used  -  the  principal  complications  come  in  the  analo- 
gous treatment  of  the  jump  conditions . ) 

As  is  known  [6],  the  simplest  explicit  difference  schemes 

for  the  Eulerlan  equations  (3),  i^) ,    (5)   are  either  unstable  or 

Inaccurate,  because  of  the  difficulty  of  treating  the  terms  containing 

the  operator  v-V   (which  here  appears  as   (v  -  x)'V,   because  of  the 

similarity  transformation).   The  MACH  code  uses  centered  space 

differences:  that  is,  a  quantity  like  —  is  approximated  at  x,  y 

^x 
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by  subtracting  values  of   f   at  x  +  Ax  and  x  -  Ax  and 
dividing  by  2Ax.   To  make  the  equations  stable,  two  sets  of 
calculations  are  made  for  each  time  step.   In  the  first  set, 
the  old  values  are  used  in  the   (v  -  x)»V  terms  to  give 
tentative  advanced  values,  denoted  by  u,  v,  p,  p,   at  the  new 
time;  in  the  second  set,  these  tentative  values  are  used  in  the 
(v  -  x)«V   terms  to  give  final  advanced  values   u,  v,  p,  p.   The 
error  is  (5'(Ax  )  +  (5'(Ay  )  +(5'(At)  (but  the  coefficients  of  the 
Cf(AT)   terms  vanish  for  a  steady  state). 
The  abbreviation 


fjj^  =  f(jAx,kAy,nAT) 


is  used  for  any  function   f(x,y,T);   n   can  be  an  integer  or 
half  an  odd  integer.   The  equations  for  the  tentative  values  are 
11  11 

J^  J^      +    (u.    2    _    jAx)    ^i^i^ i^^^ 

At  ^  2Ax 


1  1 

1  ""-2  ^-2 

^    """2         ,,    ,    ^,1    k+1    -   ^j    k-1 
(v.,       -   kAy)   -^ ^ 

3^  2Ay 


^1^    ZkjLA^lrUi  =  0,  (7) 

,1 

a  similar  equation  for  v.,   ,  i°i 


the  equation 
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Pjk  -  P.Tk 
At 


ji+  crj  -  jAx)  '"^^  "^  "  'J-^  i^ 


n- 


Jk 


2AX 


+  (v 


,1 


n 


n 


Jk 


2Ay 


+  P 


n 


Jk 


n- 


n- 


u 


u. 


,1+1  k    j-1  k 
2Ax 


+ 


,1        ,1 
^1  k+1  -  \i  k-1 


2Ay 


=  0  ,    (9) 


and 


'^.n+l    n        1 


n 


n 


At 


-  jAx) 


\1+1  k   ^j-1  k 
2Ax 


n- 


+  i'^^^^    -   kAy) 


P.1  k+1     Pj  k-1 

2Ay 


7P 


n 


+ 


Jk 


n 
'jk 


~n+l 


n 


At 


+  .  .  . 


0 


(10) 


The  terms  represented  by  the  dots  in  (lO)  are  identical  with  the 
second  and  third  terms  of  (9). 

The  equations  for  the  final  values  are  the  same  as  (7)-(lO) 

1 


,1 


n- 


except  that   l)   u.^^   ,   v.^   ,   pj^ 


p^J^  ,   and  p'^J^   are  replaced 


n- 


1 


,1 
n-l-2 


by  the  final  values  u.,  ^  ,   v  ,  ^  ,   p'^J^  and   p^J"';   in  the  first 

jK         JK         J  K  JK      -|^  2^ 


•^-2 


n- 


2 


term  of  (?),  (8),  (9),  and  (lO) ,  respectively;   2)   u  .^   ,   u  ^^ 


k 


11  - 


v^  ,      etc.  are  replaced  by  u    ,   u  ^^  ^    ,      v    ,      etc. 

,1 
In  the  second  and  third  terms  of  (7)  and  (8);   5)   u .,   ,      etc. 

,1 
are  replaced  by  u ..,    etc.  throughout  (9)  and  (lO);   4)  9 -y.   > 

p.,   etc.  are  replaced  by  p.,   ,      p..    In  {9)  and  (lO)   except 
in  the  first  terms. 

From  an  analysis  of  several  simplified  versions  of  these 
equations,  it  seems  rather  certain  that  the  scheme  is  stable  if 


^   P   2AX  P   2Ay 

throughout  the  flow;  a  proof  is  still  lacking. 

4 .  Representation  of  curved  boundaries. 

The  reflected  shock,  the  Mach  shock,  and  the  slip  surface 
are  represented  by  strings  of  points  spaced  a  little  closer  than 
the  net  points,  as  in  Figure  2.   Because  of  the  second  order 
accuracy  of  all  formulas,  this  is  equivalent  to  replacing  the  curve 
by  a  sequence  of  smoothly  joined  parabolic  arcs.   The  coordinates 
of  the  £ —  point  on  the  shock  are  called  x„      ■.,    jn      ,  ,   or 
simply  X.,  y . .   For  this  purpose  the  two  shocks  are  treated 
together  as  a  single  curve  running  counter-clockwise  from  the 
point  of  contact  with  the  wedge,  at  the  right,  to  the  point  of 
contact  with  the  floor  of  the  shock  tube,  at  the  left.   [This 
curve  has  a  corner  at  the  triple  point,  and  special  formulas  are 
used  there,  based  on  the  assumption  of  finite  curvature  of  both 
curves  at  the  corner.   It  is  furthermore  assumed  that  the  left 
terminus  of  the  curve  is  to  the  left  of  the  edge  of  the  wedge  - 
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this  Is  true  unless  the  incident  shock  Is  so  strong  as  to  make 

the  flow  behind  it  supersonic.  In  which  case  there  may  be  an  attached 

shock  at  the  wedge.]   The  representation  is  supplemented  by  giving 

the  components   >^x,i,sh,  ^y,i,sh'   °^  simply  7,^,    A  ,   of  the 
normal  unit  vector  at  each  point,  calculated  by  fitting  a  parabola 
to  the  £ —  ,   (i  +  l) —  ,   and  {£    -   l] —  points  with  its  vertex 
at  the  £ —  point.    As  the  curve  moves,  the  equations  of  fluid 
dynamics  determine  the  component  of  velocity  of  the  point  x„,  y^ 
in  the  direction  of  A   ,  and  we  are  free  to  assign  the  tangential 
component  arbitrarily;  we  do  it  so  as  to  keep  the  points  more-or- 
less  uniformly  spaced  along  the  curve.   Furthermore,  the  points 
have  initially  a  uniform  spacing  As   in  arc  length,  and  if, 
during  the  calculation,  the  spacing  anywhere  exceeds  ^  As  ,   a 
new  point  is  interpolated  and  the  others  are  renumbered.   (No 
provision  has  been  made  for  deleting  points  if  the  spacing  decreases.) 
The  interpolation  is  only  linear,  but  the  necessity  for  it  should 
cease,  any  way,  when  the  flow  is  nearly  stationary. 

The  representation  of  the  slip  surface  is  similar  to  that 
of  the  shock. 

The  provision  for  interpolating  new  points  on  the  curves 
was  added  to  the  code  because  the  slip  line  increases  considerably 
in  length  during  the  early  cycles  of  calculation. 

At  each  point  of  the  shock  or  slip-line,  the  calculation 

1    •  1 
n-i.   n-2- 

carries  along  values  of  various  quantities  like   u  ^,  v   , 

n+—   n+i 
n   n     p     p   ^^n  '^n     /       / 
p,p,u    ,v   }    P    >    p    }    dx/dT,  dy/dT  ,  A«v,   etc.   Two  complete 

sets  of  such  quantities  are  carried  on  the  slip  line  -  one  for  each 

side;  in  fact  the  slip  line  is  treated  as  two  separate  (but  coincident) 

curves  having  the  same  values  of  x,  y,  dx/dr,  dy/dT   at  its  points. 
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The  vector  A   Is  always  taken  as  the  outward  normal  from 
the  region  In  which  the  quantities  apply.   That  Is,  the  normal 
to  the  shock  Is  outward,  the  lower  side  of  the  slip  line  carries 
an  upward  pointing  normal,  and  the  upper  side  a  downward  pointing 
one, 
5.  Evaluation  of  partial  derivatives. 

For  carrying  the  solution  forward  In  time  at  any  point, 
whether  boundary  point  or  net  point,  we  need  values  of  the  partial 
derivatives  of  u,  v,  p,  p   with  respect  to  x  and  y.   At  regular 
net  points,  these  are  approximated  by  centered  difference  quotients, 
as  described  In  section  3-   At  special  net  points  and  at  boundary 
points.  It  must  be  done  differently. 

Associated  with  each  special  net  point,  like   A  In  Figure  3j 
is  a  cluster  of  8-10  points,  some  of  them  net  points  and  some 
boundary  points;  the  cluster  associated  with  A   in  the  figure  Is 
encircled  by  the  solid  closed  curve.   To  evaluate  the  partial 
derivatives  of  one  of  the  functions,  say  u(x,y),  in  the  neighborhood 
of  point   A,   the  values  of  u(x,y)  at  the  cluster  points  are  fitted 
with  a  quadratic  function  of  x   and  y,   and  the  partial  derivatives 
of  this  function  are  evaluated  at  each  of  the  points  indicated  by 
solid  circles.   The  adjacent  clusters  overlap  -  two  others  are 
encircled  by  dashed  curves  in  Figure  3-   Each  of  the  other  points 
of  a  cluster  (the  ones  indicated  by  open  circles  in  the  cluster 
about   A   in  the  figure)  receives  its  values  of   ^u/Sx  and  Su/5y 
either  from  another  cluster  in  which  it  is  more  centrally  located 
or  from  the  normal  divided-difference  procedure  used  for  regular 
net  points  -  points   B  and  C   in  Figure  3  a-ve   regular  net  points. 

The  fitting  is  done  by  the  following  least-squares  procedure. 
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If  X  ,  Y   are  the  coordinates  of  the  m —  point  of  a  cluster* 
m   m 

and  u   is  the  value  of  u(x,y}   there,  if  the  quadratic  function 
m 

2 
that  is  to  do  the  fitting  is  written  as   a^  +  a-j^x  +  a^Y   +  a^x  + 

QhXy  +  QcY  ,   and  if  we  call 


w  =  (w  )    - 
rlxa  qm 


m 


y. 


m 


m 


X  y 


V 


y. 


m 


') 


,        a  = 


O 


a-, 


v!v 


,(12) 


the  problem  is  to  minimize  the  quadratic  form 


=  y      I     \{'w    »CL   -   u)' 

^ (m)  ^jn  ^    m 


[13) 


This  has  to  be  done  for  each  of  eight  functions,  namely  the 
provisional  and  final  values  of  u,  v,  p,  p.   Much  of  the  work 
involved  is  the  same  for  each  of  the  functions;  hence,  it  is 
more  efficient  to  write  the  solution  as 


%-^ 


u 


/   \   (J      u 

-(m;   pm  m 


(1^) 


and  solve  for  the  coefficients   c    Instead  of  directly  for  the 
a  .   For  a  given  cluster  at  a  given  cycle  of  the  calculation, 


The  meaning  of  the  subscript  is  here  different  from  that  of 
section  4. 
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the  coefficients   c    can  then  be  used  eight  times  over.   The 
equation  for  them  is 


^(r.)    "Inn  C^m  =  ^nm'  (^Sl 


-(p)  "'qp  "pm  -  ^qm' 
or,  in  matrix  notation, 

MC  =  W,  (16) 

where  M  is  the  6x  6  matrix  with  elements 


m 


■qp  =   ^(r)  V  V^  ^^^^ 

cluster 


and   C  and  W  are  matrices  of   6   rows  and   8-10   columns. 

To  solve  (16)  I  might  have  used  standard  matrix  routines. 

I  might  have  used  a  linear  system  solver  repeatedly,  taking  (16) 

column-by-column;  this  is  inefficient,  because  each  of  the  8-10 

equation  systems  has  the  same  matrix  M  and  much  of  the  work 

would  have  to  be  repeated.   I  might  have  used  a  matrix  inverter 

to  give   M~   and  applied  this  to  each  row  of   W  in  turn;  this 

is  also  inefficient  because  first,  the  complete  inverse  of  M  is 

not  needed,  and  second,  the  solution  is  not  completely  needed  - 

the  constant  coefficient   a   is  not  needed  (since  we  only  use  the 

o 

partial  derivatives  of   u(x,y));  hence  the  back  substitution  may 
be  stopped  before  completion.   I  therefore  coded  a  special  routine 
to  solve  (16)  for  the  bottom  five  rows  of  the  matrix   C. 

If  the  points  of  a  cluster  should  happen  to  be  lined  up  on 
a  conic  section,  the  coefficients  of  the  quadratic  fit  would  be 
either  indeterminate  or  infinite  (they  are  infinite  unless  the  given 
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values  of  u  fit  a  quadratic  function  exactly) .   If  the  points 
lie  close  to  a  conic  section,  the  values  of  the  partial  derivatives 
obtained  by  the  fit  are  expected  to  be  Inaccurate.   The  procedure 
that  was  provided  in  the  original  code  to  take  care  of  this  con- 
tingency was  to  minimize  the  quadratic  form 

P  +  a  Ba,  (18) 

instead  of  F  alone  (eq.(l3))j   where  B  is  a  positive  semi- 
definite  matrix  whose  elements  are  supposed  to  be  small  in  some 
sense,  but  which  keeps  the  eigenvalues  of  the-  matrix  M  +  B  of 
the  form  (18)  bounded  away  from  zero.   The  choice  of  B  made  for 
the  MACH  code  will  be  given  below. 

Minimization  of  (18)  with  respect  to   a  gives 

(M  4-  B]a  =  Z](^)  w^  u^  (19) 

and  the  equation  corresponding  to  (16)  is 

(M  +  B)C  =  W.  (20) 

Hence,  the  fixed  matrix  B  is  simply  added  to  M,   given  by 
(17),  before  solving  for  C. 

Clearly,   B   should  be  so  chosen  that  the  form  a  Ba  is 
invariant  under  rotations  of  the  coordinate  system,  i.e.,  under 
the  transformations  of   (a  ,  a-,  ,  .  .  .  ,at- )   Induced  by  rotations  in 
X,  y.   A  complete  set  of  quadratic  invariants  is 
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2 


2  ^   2 

a-,   +  ttp  J 

2  ^  1   2  _^   2 


(21) 


The  form  should  also  be  invariant  under  translations  of  the  coordinate 
system,  and  this  requirement  rules  out  the  first  two  of  the  above 
expressions.   For  simplicity,   a'Ba  was  taken  to  be  proportional 
to  the  last  expression;  that  is,   B   is  the  diagonal  matrix 


^0 


B  =  K 


(o> 


0 


1 

2 


(0) 


V 


1 


(22) 


y 


where   ;<c   is  a  small  constant  (its  actual  choice  will  be  given  below). 

The  idea  behind  the  introduction  of  B   is  this:  normally 
M  +  B  %  M,   because   i<   is  small.   Specifically,  if  the  points  of 
a  cluster  are  not  close  to  a  conic  section,  the  eigenvalues  of  M 
are  not  close  to  zero,  and  hence  the  inverse  of  M  +  B   is  nearly 
the  same  as  that  of  M.   If  the  points  of  a  cluster  lie  on  a  conic 
section,  there  is  a  non-trivial  quadratic  function  ^(x,y)   that 
vanishes  at  all  points  of  the  cluster  (and  on  the  entire  conic 

section) .   If  the  points  lie  close  to  the  conic  section,  then  ^(x,y) 

2 
will  generally  be  present  in  the  minimizing  function   a  +•  •  •+  a^-y 
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■i  very  large  coefficient;  this  hardly  affects  the  value  of 
the  minimizing  function  at  any  cluster  point  but  it  can  affect  the 
values  of  the  derivatives  disastrously.   However,  the  presence  of 
Tpix.y)      will  give  a  large  value  to   a«Ba,   if  its  coefficient 
is  large:  hence,  the  presence  of  the  term  a«Ba  in  the  form 
to  be  minimized  tends  to  prevent  the  minimizing  function  from 
acquiring  an  unnaturally  large  curvature  in  its  effort  to  fit  a 
set  of  values  at  points  lying  almost  on  a  conic  section.   In  the 
limit  of  large  k,      tending  to  oo,  a^,  au,    a^     would  be  forced  to 
be  zero,  and  the  fitting  function  would  be  the  best  linear,  rather 
than  the  best  quadratic  function  of  x,  y. 

There  would  still  be  trouble  if  the  cluster  points  were  nearly 
lined  up  on  a  single  straight  line  (because  B   is  only  semi-definite), 
but  the  method  of  choosing  the  cluster  points  is  such  as  to  prevent 
this  . 

Assuming  that  upon  refinement  of  the  point  net  Ax  and  As 
decrease  at  the  same  rate,  the  size  scale  of  a  cluster  is  proportional 
to  Ax.   If  we  assume  that  for  similar  clusters  of  different  sizes, 
the  fractional  change  in  an  element  of  M  caused  by  addition  of  B 
should  be  the  same,  we  can  write 

K  =  «^(Ax)^,  (23) 

where   «:   is  a  dimensionless  constant,  presumably  rather  small 
compared  to  unity. 

One  thing  that  was  learned  from  running  the  problem,  and 
which  is  obvious  in  the  wisdom  of  hindsight,  is  that  the  matrix  M 
can  be  rather  ill-conditioned  even  if  the  points  of  the  cluster  do 
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not  approximate  a  conic  section.   This  Is  especially  the  case  If 

the  cluster  happens  to  lie  far  from  the  origin,  since  the  vectors 

X   =  (x  ,y  )   are  then  all  nearly  equal.   The  obvious  remedy  Is 

to  choose  for  each  cluster  an  origin  within  It;  a  convenient  choice 

of  origin  Is  the  position   (x  ,y  )   of  one  of  the  cluster  points; 

that  Is,   X   and  y   are  replaced.  In  (12),  by   (x   -  x  )  and 
^    m       "^m        ^  ■''''mo 

(y   -  y  );   the  approximation  to  u   is  then  a     +  a^ (x  -  x  ) 
■'m   '^o  '       ^^  o    1      o 

+  •••+  a(-(y  -  y  )  .   When  this  was  done,  no  more  trouble  with  111- 
condltionlng  was  encountered  except  for  clusters  lying  nearly  on 
a  conic  section.   Since   (x   -  x  )   and   (y   -  y  )   are  of  order 
Ax,   the  elements  of  the   3^5   submatrlx  at  the  lower  right  of 
M  are  of  order   (Ax)  ;   hence  the  choice   (23). 

A  second  thing  that  was  learned  from  running  the  problem  Is 
that  It  Is  surprisingly  easy  for  8-10  points  to  lie  very  close  to 
a  conic  section.   It  had  been  hoped  Initially  that   /<   could  be 
chosen  <<   1,   and  that  only  rarely  would  the  presence  of  the  matrix 
B   interfere  appreciably  with  the  second-order  accuracy  of  the 
quadratic  fitting.   But  experience  so  far  Indicates  that  it  is  not 
safe  to  take  k        much  smaller  than  x;  y^   ;   this  causes  considerable 
loss  in  accuracy. 

If  the  least-squares  quadratic  fitting  were  to  be  retained 
as  a  method  for  estimating  gradients.  It  would  seem  that  preferable 
to  the  Introduction  of  B  would  be  to  compute   det  M  for  each 
cluster  and  to  include  additional  points  in  the  cluster  whenever 
det  M   is  smaller  than  a  certain  value. 
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6 .  Shock  conditions. 

We  give  here  the  Ranklne-Hugonlot  jump  conditions  that  apply 
at  each  point  of  either  shock.   We  give  them  first  In  terms  of  the 
original  Eulerian  variables   X,   V  and  convert  them  to  the  similarity 
variables  later.   In  this  section,  the  symbols  V  p,  p   denote  values 
at  some  point  immediately  behind  the  shock  (inside,  in  Figure  l); 
V  ,  D  ,  p   denote  values  ahead  -  there  will  be  one  set  of  values 
ahead  of  the  Mach  shock  and  another  ahead  of  the  secondary  shock; 
in  particular,   V  =  0  ahead  of  the  Mach  shock.   The  velocity  of 
the  shock  itself  is  denoted  by  W.   Only  the  component   A'W  is- 
meaningful,  and  A*'W  may  be  chosen  at  will   (A*   is  a  unit  vector 
tangent  to  the  shock;   A*  I  a),   but  in  any  case  if  X  is  a  point 
on  the  shock  at  time   t,   then  X  +  W  At   is  a  point  on  the  shock 
at  time   t  +  At. 

The  jump  conditions  are: 


p^  A«(W  -  V^)  =  pA»(W  -  V) 


(24) 


2,*.(V  -  V  )  =  0  ,  (25) 

and 

^  =  'P  '  Po  ,  (26) 

Po    0p^  -  p 

where   0:  =  (7  +  !)/(/  -  l)- 

These  are  four  equations   ((24)   is  two  equations)   for  five 
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unknowns;   p,  p,  A*W,   and  the  two  components  of  V.   To  fix  the 
unknowns  J  an  additional  equation  Is  needed,  which  evidently  should 
specify  how  the  shock  Is  Influenced  by  the  flow  "behind  It.   For 
this  purpose  one  may  take  one  of  the  partial  differential  equations 
(3)  J  (^)j  or  (5)^  oi^j  more  generally,  a  linear  combination  of  them. 
If  the  time  derivatives  In  this  combination  are  approximated  by 
forward  difference  quotients  and  all  quantities  are  known  at  the 
earlier  time,  the  equation  then  gives  the  value  of  some  linear 
combination  of  the  above  five  unknowns  at  the  shock. 

7  .  Characteristics . 

To  decide  what  linear  combination  of  ij>) ,    (4),  and  (5)   Is 
most  appropriate  for  supplementing  the  jump  conditions  (2^),  (25), 
and  (26)  on  the  shock,  a  heuristic  appeal  Is  made  to  the  theory 
of  characteristics.   A  linear  combination  (whose  coefficients  may 
depend  on  X,  Y,  t,  U,  V,  p,  p)   of  the  partial  differential 
equations  Is  said  to  be  In  normal  form  or  characteristic  form  If 
the  directions.  In  the   X,  Y,  t   space.  In  which  U,  V,  p,  p   are 
differentiated  all  lie  in  a  single  plane.   The  possible  orientations 
of  such  a  plane,  at  a  given  point   X,  Y,  t,   can  be  obtained  by  a 
straightforward  calculation;  they  will  be  described  below.   In  a 
certain  sense,  effects  are  propagated  along  or  information  flows 
along  these  planes,  and  since  we  are  interested  in  supplementing 
(24),  (25)^  and  (26)  with  information  that  has  arrived  at  a  point 
on  the  shock  from  the  fluid  in  the  region  behind,  we  should  choose 
this  plane,  if  possible,  so  that  the  past  portion  of  it  lies  (at 
least  locally)  entirely  in  the  region  behind  the  shock.   This,  it 
turns  out,  uniquely  fixes  the  linear  combination  of  (3),  (4),  and 
(5)  that  is  to  be  used.   The  heuristic  concept  of  flow  of  information 
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suggests  that  this  linear  combination  will  couple  the  difference 
equations  for  the  smooth  part  of  the  flow  behind  to  the  jump 
conditions  at  the  shock  as  closely  as  possible;  It  thus  suggests 
that  with  any  other  linear  combination,  the  numerical  equations 
may  either  be  unstable  or  require  a  reduced  value  of  At  for 
stability. 

Equation  (5)  -  constancy  of  entropy  along  a  particle  path  - 
Is  already  in  characteristic  form;  in  this  case,  not  only  are  the 
directions  in  which  C  and  p   are  differentiated  coplanar,  they 
are  actually  the  same  direction,  namely  that  of  particle  paths  in 
the  X,  Y,  t  space.   Clearly,  this  equation  is  unsatisfactory  for 
supplementing  the  jump  conditions  at  the  shock;  since  the  shock 
overtakes  the  particles,  this  equation  yields  information  on  the 
entropy  just  ahead  of,  not  just  behind  the  shock. 

Each  component  of  equation  (3)  is  also  in  characteristic  form, 
and,  more  generally,  if  n   is  any  vector  in  the  X,  Y  plane,  the 
scalar  product  of  [i     with  equation  (5),  namely 

[I'D  V  +  -\  11. Vp  =  0,  (27) 

pS 

is  in  characteristic  form,  because  each  component  of   v  Is  differ- 
entlated  along  the  particle  path  and  p   is  differentiated  in  the 
direction  of  i_i.   In  this  case,  then,  the  plane  referred  to  in  the 
definition  of  characteristic  form  is  any  plane  tangent  to  the  particle 
path.   Plow  of  information  along  these  planes  is  also  unsatisfactory, 
because  the  past  portion  of  any  of  them  lies  mostly  ahead  of  the 
shock.   [However,  equation  (27)  will  be  used  for  slip  surface  fitting, 
described  in  section  10,  because  in  that  case  particle  paths  adhere  to 
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the  slip  line;  therefore  if  [i      is  tangent  to  the  slip  surface, 
the  plane  in  question  lies  (locally)  in  the  region  (more  precisely, 
on  the  boundary  of  the  region)  in  which  the  point  X,  Y,  t   is 
located.  ] 

Lastly,  if  i^   is  again  any  vector  in  the  X,  Y  plane,  the 
equation 


Pc(m,D  +  cV).v  +  (D  +  CM.»V)p  =  0,  (28) 


where   c   =  *  V/P/p  =  local  sound  speed,  is  a  linear  combination 
of  i^) ,    (^),  and  (5)  and  is  in  characteristic  form.   The  plane  in 
which  the  directions  of  differentiation  lie  is  now  inclined  to  the 
particle  path,  and  is  tangent  to  the  sonic  cone.   In  Figure  k, 
planes  corresponding  to   t  =  t   and   t  =  t  +  At   are  shown 
horizontally;  Q  Q'   is  a  particle  path;  the  circle  in  the  upper 
plane  about  center   Q'   has  radius   cAt,   and  the  cone  through 
this  circle  and  with  Q  as  apex  is  the  sonic  cone;  the  plane 
(j    ,      shown  in  the  figure,  is  tangent  to  the  sonic  cone.   We  must 
choose   0  so  that  its  intersection  with  the  upper  plane  is  tangent 
to  the  shock,  say  along   AB.   Then,  since  a  shock  is  subsonic  with 
respect  to  the  fluid  behind  it,  the  backward  part  of   (P  (the  part 
shown)  lies  wholly  in  the  shocked  fluid,  and  information  flowing 
upward  along  (j   overtakes  the  shock.   Therefore  yi        must  be  set 
equal  to   A,   the  shock  normal,  and  then  (28),  when  expressed  in 
difference  form,  is  to  be  used  to  supplement  the  jump  conditions. 
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8.  Shock  fitting. 

The  position  of  the  shock  and  the  quantities   v,  p,  p   at 
points  immediately  behind  it  are  now  to  be  determined,  when  all 
quantities  are  known  at  the  preceeding  cycle  of  the  calculation. 
Expressing  (24)  and  (25)  in  the  similarity  variables  (l),  we  have 


V  -  Sv,  V  =  vSv  , 


'cix 


W 


=  Sfe  +  ?.)   .  (29) 


A1 


Here  x  =  x(t)   is  one  of  the  points  x»  ^,   referred  to  in 
section  h   for  following  the  motion  of  the  shock;  in  the  present 
section,  the  subscripts  are  omitted;   v,  p,   and  p   denote 
values  of  fluid  velocity,  pressure  and  density  immediately  behind 
the  shock  at  the  point  x;    x"^,  v^,  p  ,  p   denote  their  values 
at  time  t  =  nAT . 

The  values  of  the  flow  quantities  ahead  of  the  shock  are 
denoted  by  v  ,  p  ,  p  .   Here  two  cases  must  be  distinguished:  the 
Mach  shock  and  the  reflected  shock.   The  code  distinguishes  which 
case  a  given  point  x  belongs  to  on  geometrical  grounds:  by 
suitable  choice  of  units,  the  primary  shock  is  taken  to  be  at 
X  =  1   for  values  of  y  above  the  triple  point;  therefore,  the 
Intersection  of  the  complete  shock  curve  with  the  straight  line 
X  =  1   separates  the  shock  curve  into  Mach  shock  and  reflected  shock, 
p  ,  p,  denote  pressure  and  density  ahead  of  the  primary  shock  and 
are  given.   Pp,  pp   denote  values  behind  the  secondary  shock  (above 
the  triple  point);   Pp   is  given  -  the  ratio  Pp/p-,   characterizes 
the  strength  of  the  primary  shock.   We  have: 


-  25  - 


(l)  Ahead  of  Mach  shock, 


Po  =  Pi' 


Po  =  Pi 


} 


u   =  0  , 
o     ' 


V  =  0  ; 

o     ' 


(30) 


(2)  Ahead  of  reflected  shock. 


Po  =  P2 


^P2  +  Pi  ,.        7+1) 

Po  =  P2  =  Pi : ^^  =  T^^^' 

0Pl  +  P2 


(31) 


U    =  —    /(p^  -  Pn  )  ( -  )    ) 

O    ^      \/'^^  ^1   p^     P2 


V   =  0. 

o 


(32: 


At  shock  points   the  calculation  carries  along  the  values 
of  p   and  p   for  Integer  n   and  of  u  and   v  for  half-odd- 
integer  n,   just  as  at  net  points;  in  addition,  it  carries  along 
values  of  p   for  half -odd-integer  n  and  of   A«v  for  integer  n, 
Values  of  x   and  y  are  carried  for  integer  n,   and  values  of 
dx/dT   and  dy/dT    for  half -odd-integer  n.   The  shock-fitting 
formulas  take  four  slightly  different  forms  according  to  which  of 
the  following  sets  of  values  are  to  he  obtained: 


-  26  - 


1.  tentative  values  at  n  +  ^  , 

2.  tentative    "    "  n  +  1  , 

3.  final 
i|.  final       "    "  n  +  1  . 


"    "   n  +  I  , 


Only  the  first  forai  will  be  described  here,  as  the  others  are 
similar.   The  changes  for  going  from  tentative  to  final  values  are 
of  the  same  kind  as  described  at  the  end  of  section  3  for  regular 
net  points. 

To  calculate  preliminary  values  at  time  (n  +  p-)AT,  when 
all  quantities  are  known  up  to  time  nAx,  the  finite-difference 
equivalent  of  (28)  Is  written,  with  p.    taken  as  A  ,   and  with 


D   written  out  as  -5—  +  (v  -  x  -  dx/dx )  •  V  [-=—  denotes  time 


•dT 


derivative  of  a  function  of  x  and  t  after  x  has  been  set 
equal  to  x(t),   the  coordinates  of  the  point  on  the  shock].   The 
equation  is 


Sn  n 
c  p 


1 

,  n    2 
A  .  V 


,n  ""^       ,n  ,^^-2        .  ^  1         .      n-i   .  n  ^^ 
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(35) 
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n-l-2    n-g- 
P     -  P 

AT 


dx   2 


/  n-7s-         dx   ^  "N 


1   n^n  „  n    „ 
+  p^  c  A  •  Vp   =0; 


it  is  solved  to  give  the  value  of  the  quantity 
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,1  ,1  .1 


■a     2      def.  2    ,    C  n  Yif^n  2\  /^),\ 

P  =P  +Ocp(A.v),  (34) 

In  which     p  and     A'V     are  unknown.       (All   other  quantities 


r\j      r\j 


In  (33)  are  known.) 

Equation  (3^)  and  the  four  Ranklne-Hugonlot  Jump  conditions 

are  solved  simultaneously  (by  a  Newton-Raphson  procedure)  for   -, 

,1      ,1         ,1  .1  ^   n^ 

n+2     n+g-        n-F^-         n+g-  dx    2 

the  unknowns   p    ,   p    ,   (A»v)    ,   (A*»v)    ,   and  i'K'-P^) 
Here,   A   stands  for  A  ,   and  A*   stands  for  the  tangent  unit 
vector  at  time   nAT .   In  the  present  notation  (see  equation  (29)), 
the  jump  conditions  are 


1  /  ,1 

n-^-^  /   n- 


Si-(v  "  -  vj  =y(p  ^  -p„)(i ^-j)        ,  (36) 

^o    n-+2- 
P   ■ 

n^ 
A^'(v  2  _    )  =  0   ,  (37) 


and 


,1 
n+2      0p     -  pQ 

P   =  Po T 

n+2 
0Po  -  P 


(38) 


Subsequent  to  the  original  coding,  equation  (3^)  was  replaced 
by  the  equation 

,1      ,1  n-f^ 

„   2      2  ,  C  n  n  A  -v    ) 
P    =  p    +  Oc  p   '%.  'X, 

(34') 
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where  \i.      Is  a  positive  constant.   The  added  term  vanishes  as 
0((As)  )   on  refinement  of  the  net^  hence  does  not  change  the  formal 
order  of  accuracy;  it  helps  to  stabilize  the  shock,  because,  when 
written  as 

(As)2   -I±l 2ll 2lirl  .  ^, 

^  (As)2 


it  represents  a  surface  tension  (compare  with  equation  (42)  in 

Hi 


o 

section  10)  with  coefficient  u.-,  (As)  . 


,1 
n-Fp- 

In  order  to  solve  (33)  for   P    ,   the  spatial  derivatives 

are  obtained  as  described  in  section  5' 

The  calculation  of  the  final  values  at  time   (n  +  p-)At  differs 

1  ,1 

from  the  above  in  that  v     is  replaced  by  v     throughout  (33) 

except  for  the  first  term  in  the  square  bracket.   The  calculations 

of  the  preliminary  and  advanced  values  at  time   (n  +  1)At   are  similar. 

9 .  Motion  of  the  shock. 

As  soon  as  the  final  values  at  time   (n  +  p-)AT   have  been 

obtained  on  the  shock,  the  shock  is  moved  to  its  new  position. 

Convenient  (small)  values  are  assigned  to  the  tangential  motion 

,1 
A*»(dx/dT)  "^      ,  (59) 


f\/ 


so  adjusted  as  to  tend  to  eliminate  irregularities  of  spacing  of 
the  points  along  the  shock,  and  then  the  new  position  is  found  from 
the  equation 

dx  n-h=r      ,  n 

XH  +  At(;^)    2  =  x''^^   .  (40) 
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[Actually  a  slight  refinement  of  this  procedure,  perhaps  of  doubtful 

value,  was  made  in  the  code.   Namely,  In  an  effort  to  center  (35) 

more  accurately  at   (n  +  75-)^   x   was  replaced,  in  the  left  member, 

^    ^     dx^^  n+i 

by   X  +  7=r  Ai  {-^)         ,      before  solving  for   (dx/dx)    .] 


After  all  points  x.  \i    =   0,1,2,...)   of  the  shock  have  be 


en 


moved,  in  accordance  with  (40),  a  new  normal  vector  A   is  computed 

at  each  point  as  the  normal  to  a  parabola  passing  through  x„      , 

'^Z  -1 

x„,      and  x„,T   and  with  its  vertex  at  x.. 

The  boundary  conditions  are  that  the  bottom  point  of  the  Mach 
shock  is  constrained  to  move  along  the  wedge  and  to  have  its  normal 
vector  pointing  upward  along  the  wedge  and  that  the  bottom  point 
of  the  reflected  shock  is  constrained  to  move  along  the  floor  of 
the  shock  tube  and  to  have  its  normal  vector  pointing  horizontally 
to  the  left. 

10.  Slip-line  fitting. 

The  fitting  procedure  for  the  slip  line  may  be  contrasted 
with  shock  fitting,  as  follows:  l)  there  are   9   unknowns  instead 
of   5^   because  the  flow  quantities  are  unknown  on  both  sides; 

2)  there  are   3   jump  conditions  instead  of   4,   because  the 
tangential  component  of   v   is  no  longer  required  to  be  continuous; 

3)  on  each  side,  there  are   3  acceptable  linear  combinations  of 
the   4   partial  differential  equations  to  supplement  the  jump 
conditions.  Instead  of  just   1. 

Flow  quantities  on  the  upper  and  lower  side  of  the  slip  line 
will  be  denoted  by  superscript   "+"  and   "- "  ,   respectively. 
The  unknowns  are  the  eight  quantities   u—  ,  v—  ,  p—  ,  p—  and 
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the  normal  component  of  the  velocity  of  the  slip  line  Itself. 

The  kinematic  jump  condition  is  that  the  normal  component 
of  V  is  the  same  on  both  sides  and  is  equal  to  the  normal 
component  of  the  velocity  of  the  slip  line.   By  {29),    this  gives 
the  two  equations 

A.v  =  A.v  =A'(-^+x),  (4l) 

where  x  =  x(t)   is  one  of  the  slip-line  points.   In  (4l),   A 

+        -  -      + 

can  be  either  A   or  A  ,   since  A   =  -A  . 

Because  of  Helmholtz  instability,  discussed  below,  a  surface 

tension  has  been  attributed  to  the  slip  line.   [This  is  unphysical 

(at  least  if  the  fluid  is  an  ideal  gas),  but  the  amount  of  surface 

tension  required  is  quite  small,  and  the  hope  was  that  after  the 

flow  has  settled  down  to  a  steady  state,  the  surface  tension  can 

be  removed  (perhaps  gradually)  for  the  last  several  cycles  of 

calculation.   In  the  observed  Mach  reflection,  the  slip  line  is 

very  nearly  straight  and  Helmholtz  instability  should  be  slow  to 

develop.]   Consequently  there  is  a  pressure  jump  proportional  to 

the  local  curvature  of  the  slip  line,  the  pressure  being  higher 

on  the  inside  of  the  curve.   Therefore, 

,  A"'  (x.    -  2x.  -  X,   ) 

P     P    I.,  ^^^^2  '  ^^-^ 

whei^e  [i        is  the  surface  tension  coefficient  and  As   is  the  spacing 
of  points  along  the  slip-line.   The  choice  of  \i        will  be  discussed 
in  the  next  section. 

Equations  (4l)  and  (42)  must  be  supplemented  by  Information 
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from  the  partial  differential  equations  on  either  side,  to  fix 
the   9   unknowns;   6   further  relations  are  needed.   Two  of  these 
are  obtained  from  characteristic  forms  In  which  differentiations 
are  In  the  planes  tangent  to  the  sonic  cone  that  slant  toward 
the  slip-line,  one  from  each  side.   As  In  the  case  of  the  shock, 
these  give  values  of  the  quantities 

P    =   p   +  0(pc)   A  .  V  ; 

and  (43) 

p-  def.   p-  ^  5(pc}-  A".y", 


O    '^       '\j 


,+ 


(pc)—  being  values  corresponding  to  a  slightly  earlier  time 


o 


(earlier  by  p^  or   At)   and  hence  known. 

Two  other  pairs  of  partial  differential  equations  In 
characteristic  form  must  be  used.   One  pair  expresses  constancy 
of  entropy  along  streamlines  on  each  side: 

where  -p-  means  -^  p(x(t  )  ,y(T )  ,t  )   and  similarly  for  ^  . 

Equations  (4l),  (42),  (4^)  give  the  pressures,  the  normal 
component  of  fluid  velocity,  and  dx/dx .   The  finite-difference 
form  of  (44)  gives  the  densities  In  terms  of  the  pressures.   The 
tangential  components  of  the  fluid  velocity  are  still  needed,  and 
are  obtained  from  equation  (27;  applied  on  each  side,  where  \i 
Is  the  tangential  vector  A*;  this  equation  Is  In  normal  form, 
the  directions  of  differentiation  lying  in  a  plane  tangent  to 
the  particle  path  and  to  the  slip  line;  it  gives  the  tangential 
component  of  acceleration: 

-  32  - 


^"•[If -"  ^z-^- ^^'iz ■" ^2 r-^p  =  °-  (^5) 


•-u 


dx 

In  finite-difference  form,  (^5)  gives  a  new  value  of  the  component 

of  V  in  the  direction  A*   in  terms  of  the  old  value  of  its  com- 

ponent  in  the  same  direction  and  intermediate  values  of  the  pressure 

gradient . 

Just  as  in  the  net  point  calculation  and  in  shock  fitting, 

these  formulas  are  used  twice.   First,  old  values  of  v  are  used 

in  the   v»V  terms  to  give  preliminary  advanced  values  of  the  flow 

quantities,  and  second,  the  new  preliminary  values  of  v  are  then 

used  to  give  final  advanced  values  of  the  flow  quantities.   Treating 

the  boundaries  consistently,  in  this  respect,  with  the  net  points, 

where  the  difference  equations  are  believed  to  be  stable,  was  expected 

to  contribute  to  the  stability  of  the  overall  system. 

The  procedure  for  moving  the  slip  line,  in  each  cycle,  is  the 

dx 
same  as  that  for  moving  the  shock,  except  that  A»^  is  now  obtained 

from  (4l) . 

11.  Slip  line  fitting;  alternate  procedures  for  the  spatial  derivatives 

Two  procedures  are  available  for  evaluating  the  spatial  deriva- 
tives in  (44)  and  (45).   In  the  original  code,  they  were  obtained 
from  the  least-squares  quadratic  fit  to  function  values  at  points 
of  clusters,  as  described  in  section  5' 

However,  the  spatial  derivatives  in  (44)  and  (45)  are  in 
directions  tangent  to  the  slip  line,  and  an  alternative  procedure 
is  to  form  difference  quotients  with  respect  to  arc  length  of  values 
on  the  slip  line.   It  is  unnecessary,  and  possibly  even  unwise,  to 
Include  any  values  out  in  the  flow  away  from  the  slip  line  in  the 
difference  formula:  these  values  ought  to  cancel  out  in  (44)  and 
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(^5),  but  probably  don't,  quite,  because  of  the  approximate  nature 
of  the  difference  formula.   [These  remarks  do  not  apply  to  the 
formulas,  like  (33),  from  which   P   and   P~   are  obtained;  there, 
values  away  from  the  slip  line  are  necessarily  Involved,  because 
differentiations  are  In  a  plane  tangent  to  the  sonic  cone.] 

According  to  (4l),  the  vector  v   -  x   -   dx/dx   is  in  the 
direction  of  the  tangential  unit  vector   A*.   Therefore,  if   s 
denotes  arc  length,  we  can  write 


dx 
(v  -  X  -  :5^).V  = 


o^      -"w 


dT 


dx 
A*«(v  -  X  -  -T^) 

'\^  r^  r^  QT 


d_ 
ds 


(46) 


wherever  it  appears  in  (44)  and  (45).   [The  conventions  used  in 
the  code  are  such  that   A*   points  in  the  direction  of  Increasing 

s.] 

Differentiation  with  respect  to   s   is  done  by  taking  centered 
difference  quotients,  corrected,  to  second  order,  for  possible  unequal 
spacing  of  the  points  along  the  slip  line. 

I  now  believe  that  this  alternative  is  probably  superior  to 
the  first,  although  the  experimentation  with  the  code  is  so  far 
inconclusive.   A  variant  of  the  code  containing  this  alternative 
was  prepared,  but  tried  out  only  in  the  SUBMACH  problem  (see  section  20) 

12 .  "Topology"  routine. 

The  shocks,  the  slip  line,  and  the  surface  of  the  wedge  divide 
the  half  plane   y  ^  0   into  several  regions;  these  shift  somewhat 
from  cycle  to  cycle  because  of  the  motion  of  the  reflected  and  Mach 
shocks  and  the  slip  line.   After  each  such  shift,  it  is  necessary 
to  determine,  for  each  net  point,  which  region  the  net  point  lies 
in  and  attach  an  appropriate  region  index  to  it.   This  is  done  by 
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a  routine  that  uses  the  following  elementary  notion.   Let   C  be 

a  simple  closed  p     :i;   Let   P  be  a  point  on  it,  let  L  be 

a  half  line  extending  horizontally  from  P   to  the  right  to 

X  =  +  oo;  let  P  move  around  the  curve  once.   Then  L  sweeps  up 

and  down  over  a  certain  area,  and  a  given  net  point  A  lies  inside 

or  outside  the  polygon  C  according  to  whether  L  sweeps  over  the 

point  A  an  odd  or  even  number  of  times . 

In  the  routine,  {rl)  denotes  the  region  index  of  net  point 

(j,k)   and  ^   a  fixed  constant.   As   P   is  moved  along  one  side  of 
the  polygon  at  a  time,  certain  inequalities  tell  whether  net  point 
(j,k)   is  swept  over.   For  every  net  point  that  is  swept  over,  the 
region  index  is  modified  according  to  the  substitution 

After  the  entire  polygon   C  has  been  described,  the  region  indices 
for  points  outside   C  are  left  unaltered.   By  setting  all  region 
indices  initially  to  zero  and  using  different  values  of   ^   for 
different  polygons,  it  is  clearly  possible  to  assign,  in  this  way, 
a  unique  index  to  each  region. 

In  this  routine,  the  shocks  and  the  slip  line  have  been  repre- 
sented by  polygonal  curves  rather  than  the  sequence  of  smoothly 
joined  parabolic  arcs  referred  to  in  section  ^,  merely  to  simplify 
some  of  the  inequalities  that  have  to  be  tested.   This  clearly  makes 
no  difference  unless  a  net  point  happens  to  lie  between  the  parabolic 
arc  and  its  chord,  and  this  is  regarded  as  a  very  improbable  occurrence, 
because  the  arcs  are  very  short. 
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since  most  net  points  will  usually  have  the  same  region 
index  from  cycle  to  cycle,  one  ought  on.ly,  in  principle,   to 
have  to  examine  those  net  points  that  lie  fairly  near  the  shock 
or  slip  line  before  it  is  moved.   It  is  my  impression,  however, 
that  it  would  take  more  machine  time  to  decide  which  points  to 
examine  than  it  does  simply  to  examine  them  all. 

13.  Initial  conditions. 

I  originally  supposed  that  I  could  be  rather  cavalier  about 
the  initial  conditions,  because  only  the  asymptotic  steady  state 
(relative  to  the  similarity  variables)  is  of  interest,   I  supposed 
that  the  reflected  shock  could  be  represented  by  any  convenient 
circular  arc  meeting  straight  lines  for  the  primary  shock,  the 
Mach  shock,  and  the  slip  line  and  that  values  of  u,  v,  p,  p 
could  be  assigned  any  convenient  smooth  values  in  the  flow. 

The  early  cycles  of  the  calculation  showed,  however,  that 
unless  the  Initial  conditions  are  rather  carefully  chosen,  there 
are  rather  violent  transient  disturbances  that  make  checking  out 
of  the  method  difficult.   In  particular,  failure  to  satisfy  the 
Rankine-Hugoniot  conditions  on  the  shock  causes  either  a  rarefaction 
or  a  secondary  shock  to  proceed  back  into  the  shocked  material . 

New  set-up  routines  were  therefore  provided,  so  as  to  give 
initial  data  that  satisfy  all  jump  conditions.   The  various  jump 
conditions  would  be  in  conflict  at  the  triple  point  unless  the 
conditions  of  the  three-shock  configuration  (see  [l])  were  satisfied; 
hence,  a  routine  was  also  provided  for  generating  these  configurations 
This  routine,  which  will  be  described  in  the  next  section,  treats  the 
immediate  neighborhood  of  the  triple  point  and  assumes  that  each 
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boundary  may  be  regarded  as  plane  In  the  neighborhood  and  that 
the  flow  is  constant  between  these  planes. 

The  manner  in  which  the  three-shock  configuration  is  incorpo- 
rated in  the  initial  data,  a  somewhat  arbitrary  procedure,  can  be 
explained  with  the  aid  of  Figure  5«   Initially,  the  Mach  shock  and 
the  slip  line  are  assumed  straight  and  the  secondary  shock  is  assumed 
circular.   The  Mach  shock  is  taken  perpendicular  to  the  wall  of  the 
wedge,  to  satisfy  the  condition  that  the  flow  behind  it  must  be 
parallel  to  the  wall .   The  values  of  the  angle   p  and  of  the 
primary  shock  strength  Pp/p-,   are  supplied  to  the  three-shock  sub- 
routine, which  then  computes  the  angles   5-,   and  5^  and  the  flow 
quantities  V,  p,  p  in  each  of  the  four  regions.   V  is  the  fluid 
velocity  in  a  frame  of  reference  in  which  the  triple-point  is 
stationary.   In  the  frame  in  which  the  fluid  ahead  is  stationary, 
the  velocity  in  region  2  is  V^  -  V,  ,   or,  in  the  similarity  values, 
v^  =  (Vo  -  V-,  )/,C  (this  was  called  v  in  section  8  on  shock  fitting), 
v„   is  in  the  x  direction  and  has  magnitude  u^ .   The  circular  arc 
representing  the  secondary  shock  is  centered  at  x  =  Up,  y  ==  0,  as 
it  would  be  if  it  were  spreading  out  at  the  same  speed  in  all 
directions  through  the  uniformly  moving  fluid  in  region  2;  its 
radius  has  then  to  be  chosen  so  that  the  angle   &-,   between  the 
primary  and  secondary  shock  is  in  agreement  with  the  three-shock 
theory.   The  remaining  construction  is  a  simple  geometrical  one. 

It  is  assumed  that  p  and  p   have  initially  the  same  values 
Immediately  behind  each  shock  along  its  entire  length  as  in  the 
vicinity  of  the  triple  point.   The  Rankine-Hugoniot  jump  conditions 
then  determine  the  values  of  v  Immediately  behind  each  shock  and 
the  normal  component   A»dx/dt   of  the  velocity  of  the  shock  itself. 
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From  the  latter^  the  Initial  values  of 

dx   .  dX 

dk  "  S  dt  ~  ^ 

are  obtained  along  each  shock,  upon  setting  the  tangential 
component   A*»dx/dT   equal  to  zero. 

The  Initial  values  of   v   on  the  slip  line  and  in  the  interiors 
of  regions  3  and  4  must  be  chosen  as  smooth  functions  of  x  and  y 
that  satisfy  the  various  boundary  conditions.   v  might  have  been 
taken  constant  along  each  side  of  the  slip  line,  but  it  was  thought 
to  be  more  realistic  to  make  the  difference   v  -  v^  decrease  somewhat 
to  the  left.   This  was  done  by  taking 

V  =  v„  +  f(v  -  v^)  ,  (48) 

triple  point 

where   r   is  distance  from  the  center  of  the  circle  that  represents 
the  secondary  shock,   r  =  \/(x  -  Up )   +  y   .   Because  the  same  factor 
r/R  is  used  on  both  sides,  the  Jump  condition  is  satisfied  all  along 
the  slip  line,  since  it  was  satisfied  at  the  triple  point  by  the 
three-shock  subroutine. 

The  remaining  boundary  condition  on   v   is  that  its  normal 
component  must  vanish  on  the  floor  of  the  shock  tube   (y  =  0,  x  ■==  O) 
and  on  the  wall  of  the  wedge   (y  =  x  cot  p.,  x  >  O)  .   Smooth  initial 
values  of  u  ' and   v  were  assigned  in  the  interiors  of  regions   3 
and  4  by  the  simple  but  arbitrary  assumption  that  they  satisfy 
Laplace's  equation  subject  to  the  boundary  conditions  stated;  they 
were  obtained  at  net  points  by  a  standard  (extrapolated  Liebmann) 
relaxation  calculation,  which  is  part  of  the  set-up  routine  but 
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which  utilizes  much  of  the  machinery  of  the  main  code  for  differ- 
entiation at  regular  and  special  net  points  and  for  satisfying 
the  boundary  conditions  (those  on   the  rigid  walls  are  described  in 
section  1^,    below.) 

Pressure  p  and  density  p  were  taken  constant,  initially, 
in  each  of  the  regions  3  and  4. 

These  initial  conditions,  it  was  found,  send  the  calculation 
off  to  a  smooth  start. 

l4.  Three-shock  configurations. 

As  in  section  13  above,   V,  p,  p   with  subscript  a  =  1,2,3^4 
denote  values  of  the  flow  quantities  in  region  a,   in  a  frame  of 
reference  in  which  the  triple  point  is  stationary.   From  the  main 
code,  values  of   e(:=(7  +  l)/(7  -  l))^  P]_^  P^^  P2>  P2^  ^2'  ^^  ^ 
are  available.   Values  of  V  must  be  found  in  each  region,  of  p 
and  p    in  regions   3  and  4,   and  of  the  angle   6-,,   by  solving 
the  jump  conditions  on  the  three  shocks  and  the  slip  line.   The 
primary  shock  is  parallel  to  the  y-axis. 

The  steps  are 

1)    -  5  -^  ^ix 

2  )  5^2    ~    ^  ^  ^2x 

3)   assume  trial  value  of  V-, 


ly 


4 )  V,  ^  v„ 
ly   -^  2y 


5)   Obtain  p^,  p^  by  solving  the  equations 
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Px     Qp-^  -  p 


1 


Pi   0pl  -  P3 
6}  P3  ->  P4 


^     0P4  +  Pg 

^^   P2  0P2  +  P4  -^n 


1)   Obtain   5.   by  solving 


"^2x  °°^  ^1  "*"  "^2v  ^■^'^  ^1 


9)  v^  +y(p7T^^HirTl_](eos  p,  sm  p)  ->  y. 


10)  y^  +y.p^  -  p^)  j^-^)(cos  5,  sm  6)  ^  y^ 

11)  Test  whether  V_  ||  Vl;   that  Is,  compute 

return  to  step  3)  and  adjust  trial  value   of  V-,    until 

Y  =  0   to  desired  accuracy.   The  adjustment  was  made  according 

to  the  Newton  Iterative  procedure. 

At  step  8),  two  values  of   5-,   can  be  obtained.   I  believe, 
on  heuristic  grounds,  that  the  larger  of  the  two  will  give  an  Initial 
configuration  closer  to  the  desired  asymptotic  one,  but  that  both 
choices  of   6-,   will  lead  to  the  same  result  asymptotically  for 
large   t.   This  presumes  that  there  Is  only  one  possible  asymptotic 
configuration  (or  at  least  only  one  stable  one),  toward  which  all 
solutions  tend.   Questions  of  this  kind  can  be  Investigated  when  a 
successful  computing  method  is  available.   This  ambiguity  of  the 
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solutions  of  the  three-shock  problem  Is  discussed  by  Bleakney  and 
Taub  [l],  but  their  discussion  does  not  appear  to  Indicate  which 
value  of   5-,   we  should  take  in  our  calculation. 

15-  Rigid-wall  boundary  condition. 

The  rigid  wall,  consisting  of  the  floor  of  the  shock  tube 
and  the  surface  of  the  wedge,  might  have  been  treated  as  an  interface 
on  one  side  of  which  the  material  is  permanently  at  rest,  hence  subject 
to  the  boundary  condition 

,      dx 
A.v^  =  A»(;^  +  x)  =  0  149] 

in  place  of  (4l).   That  is,  the  wall  might  have  been  represented 
as  a  string  of  points,  at  which  values  of  u,  v,  p,  p   are  obtained 
by  a  procedure  like  that  on  one  side  of  the  slip  line.  [F^om  the 
point  of  view  of  the  fitting  procedure,  an  interface  is  the  same 
thing  as  a  slip  line.] 

Instead,  a  reflection  principle  was  used  on  each  straight  portion 
of  the  wall.   The  computational  point  net  was  so  chosen  that  the  floor 
of  the  shock  tube   (y  =  0,  x  <  O)   was  a  line  of  net  points.   Point 
(j,k)   has  coordinates   jAx,  kAy,   and  the  reflection  principle  used 
for  k  =  0,  J  -«  0  was  that  values  of  u,  v,  p,  p    can  be  assigned 
at  the  virtual  net  points  with  k  =  -1  (images  of  points  with 
k  =  +l),   according  to  the  equations 

u(x,-y)  =  u(x,y}, 

v(x,-y)  =  -v(x,y )):,,,, 

(50) 
p(x,-y)  =  p(x,y)  , 

p(x,-y)  =  p(x,y)  ; 
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with  these  assignments,  calculations  can  be  made  for  k  =  0 
just  as  for  ordinary  Interior  net  points. 

The  leftmost  net  point  on  the  line  y  =  0,   just  behind 
the  shock.  Is  exceptional,  because  Its  left  nearest  neighbor 
Is  not  In  the  same  region  with  It.   It  Is  treated  like  any  other 
special  net  point,  as  described  In  section  ^,    except  that  the 
cluster  associated  with  It  contains  the  reflection  In  the  line 
y  =  0  of  each  of  Its  points,  at  which  values  are  assigned 
according  to  (50). 

The  surface  of  the  wedge  Is  treated  according  to  the  same 
reflection  principle,  but  the  details  are  more  complicated,  because 
this  surface  Is  not  a  line  of  net  points.   Values  of  u,  v,  p,  p 
are  not  calculated  on  the  surface;  any  net  point  just  above  the  wall 
Is  treated  as  a  special  net  point,  as  described  In  section  5^  except 
that  now  the  cluster  associated  with  it  contains  reflections  In  the 
wall  of  Its  points  (see  Figure  6)  according  to  the  rules. 


net  point 


< 


X 

y 

u 

V 

p 
p 


jAx 
kAy 


u  ., 

"jk 
Pjk 
Pjk 


(51) 
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/ 


D 


X 


def 


jAx  cos  P  +  kAy  sin  p 
2D  cos  P  -  jAx 


image : 


< 


y  =  2D  sin  P  -  kAy 

def 
E   =   u.,  -cos  p  +  v.,  sin  g 
Jk     ^    jk 


u 


2E  cos  P  -  u 


(52) 


V  =  2E  sin  p  -  V 


P  =  P 


P  =  P 


Jk 
jk 


V 


In  case  such  a  cluster  contains  any  slip  line  points  or  any  shock 
points  (it  may,  if  it  is  sufficiently  close  to  the  point  of  contact 
of  the  slip  line  or  of  the  shock  with  the  wedge),  it  contains  their 
images,  too. 

The  slip  line  meets  the  wedge  at  such  a  small  angle  that 
between  them  there  is  generally  a  considerable,  thin,  V-shaped 
region  devoid  of  net  points.   The  slip  line  points  in  such  a  region 
need,  of  course,  to  be  included  in  some  clustery  hence,  special 
clusters  like  the  one  shown  in  Figure  7  are  allowed  for.   For 
purposes  of  cross-referencing,  this  cluster  is  formally  associated 
with  the  net  point   "  X  "   (which  actually  lies  in  the  wall), 
although   "  X  "   is  not  a  member  of  the  cluster. 

16.  Arrangement  of  the  calculation. 

A  summary  of  the  calculation  is  given  in  Table  I.   The  long 
list  referred  to  is  a  sequence  of  storage  blocks,  one  for  each  of 
the  cluster  of  points  described  in  section  5j  each  block  contains 
7^   items,  including  identification  of  the  members  of  the  cluster, 
the  coefficients   c   ,   and  the  partial  derivatives  of  u,  v,  p,  p 
at  the  special  net  point  with  which  the  cluster  is  associated  (the 
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partial  derivatives  at  the  other  points  of  the  cluster  are  stored 
elsewhere).   The  pictures  referred  to  are  outline  drawings  plotted 
on  the  cathode  ray  tube  of  the  704  and  photographed.   The  permanent 
clusters  referred  to  are  clusters  like  the  one  shown  at  the  right 
In  figure  6,   Each  net  point  Immediately  above  the  wedge  has  a 
permanent  cluster  permanently  associated  with  It  (since  It  Is 
economical  to  compute  these  clusters  at  the  beginning  of  the 
calculation)  even  though  the  cluster  will  not  be  used  during  those 
cycles  of  the  calculation  when  a  shock  or  slip  line  cuts  across  it; 
during  these  cycles  an  ordinary  cluster^  like  the  one  at  the  left 
in  Figure  6,    will  also  be  associated  with  this  net  point  and  will 
be  used  in  the  calculation;  it  will  be  described  in  another  block 
of  the  long  list.   Except  for  the  permanent  clusters,  the  long  list 
is  revised  every  cycle. 

The  entropy  and  Bernoulli  tests  referred  to  are  the  following. 
During  the  phase  called  output  B,  a  family  of  streamlines  Is  con- 
structed (these  are  the  paths  that  would  be  followed  In  the  similarity 
variables   x  y  by  fluid  elements  in  a  steady  flow  with  the  currently 
available  velocity  field  u,v).   Along  a  streamline,  the  following 
quantities  should  be  constant  when  a  steady  state  is  reached. 

-7 


a)   pp  ^  =  C^  (53i 

2 

^    yp 

(7  -   l)pb    J 


b)  '^ ;;  ^'  + — ^p  _^  +  f  iv  -  xids  =  c,       {54} 


The  latter  Is  the  analogue  of  Bernoulli's  law;  the  unfamiliar  terms 
come  from  the  fact  that  the  flow  is  not  stationary  in  the  original 
Eulerlan  variables   X,  Y,   but  in  the  similarity  variables  x,  y. 
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T:     nstancy  of  these  quantities  is  tested  to  give  an  Indication 
of  the  accuracy.  .  ;.  . 

17.  Results  to  date. 

As  stated  In  the  Introduction^  the  results  to  date  are  incon- 
clusive.  A  numerical  Instability,  not  yet  explained,  develops  on 
and  near  the  slip  line  after  5O-6O  cycles  of  calculation.   The  flow 
pattern  has  not  settled  down  sufficiently  by  this  time  to  permit 
quantitative  test  of  the  method.   Figure  8  shows  a  succession  of 
configurations  in  a  typical  run.   I  conjecture  that  the  instability 
can  be  eliminated  by  an  Improved  method  of  estimating  partial  deriva- 
tives, as   described  in  part  C  of  this  report.   Other  parts  of  the 
calculation  appear  to  run  satisfactorily,  but  this  still  has  to  be 
proved,  and  the  degree  of  accuracy  has  to  be  investigated.   The 
method  described  in  section  5  for  approximating  partial  derivatives 
was  tested  by  applying  it  to  analytic  functions  of  x  and  y  and 
comparing  with  the  exact  derivatives.   The  results  usually  agreed 
to  a  small  fraction  of  a  percent,  but  occasionally  the  disagreement 
reached  1-2  /o   for  certain  clusters.   I  was  not  able  to  Improve 
this  situation  appreciably  by  adjustment  of  the  constant   /<  . 

PART  B  The  SUBMACH  code. 

18.  Introduction. 

Because  of  the  instability  associated  with  the  slip,  line  in 
the  MACH  calculation,  I  decided  recently  to  study  a  simplified  problem, 
in  which  a  slip  line  is  the  main  feature,  and  for  which  an  analytic 
solution  is  available,  at  least  in  a  certain  approximation. 

We  consider  a  nearly  plane  interface  with  a  sinusoidal  corru- 
gation of  small  amplitude.   In  zeroth  approximation,  the  flow  is 
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uniform  on  each  side,  the  jump  in  fluid  velocity  at  the  Interface 
being  directed  across  the  corrugations.   To  make  conditions  similar 
to  those  in  the  MACH  problem,  the  compressibility  of  the  fluid  and 
the  surface  tension  of  the  interface  were  taken  into  account,  but 
the  jump  of  the  zero-order  density  was  ignored.   To  make  the  problem 
finite  for  the  numerical  work,  rigid  walls  were  placed  parallel  to 
the  interface  on  either  side.   This  leads  to  the  problem  depicted 
in  Figure  9-   A  frame  of  reference  is  chosen  in  which  the  zero- 
order  velocities  are  equal  in  magnitude  and  oppositely  directed. 
A  periodicity  condition  applies:  all  quantities  have  the  same  values 
at  X  =  1   as  at  x  =  0. 

An  initial  configuration  and  initial  values  of  u,  v,  p,  p 
are  given,  and  the  problem  is  to  find  the  motion  subsequently. 
Depending  on  the  amount  of  surface  tension  present,  the  corrugation 
of  the  interface  may  oscillate  or  grow  exponentially.   From  a 
linearized  analytic  calculation,  described  in  the  next  section, 
initial  data  are  obtained  that  should  lead  to  a  pure  oscillation 
or  pure  exponential  growth,  so  long  as  the  amplitude  is  small.   The 
numerical  and  analytic  solutions  are  then  automatically  compared, 
at  subsequent  cycles,  by  the  computer. 

19 .  Analytic  solution  of  the  linearized  problem. 

We  derive  formulas  for  Helmholtz  instability  as  modified  by 
l)  compressibility,  2)  surface  tension,  3)  the  presence  of  the 
rigid  walls . 

The  position  of  the  interface  is  described  by  the  equation 

y  =  g(x,t)  «  y  ,  —  «   1,  (55) 

°   ax 


-  46  - 


X 


H'V-^ 


(56) 


We  assume  Isentroplc  potential  flow; 


p  =  Ap^  ,  c  =  TPq/Po  ; 


(57) 


Vxv  =0,v=v  +V0=(^); 


V 


for   y  >•  g(x,t) 
for   y  <   g(xjt) 


The  velocity  potential  is 


0  =  0(x,y,t)  = 


0_^(x,y,t)   for  y  >  g(x,t) 


0  (x,y,t)   for  y  <     g{x,t) 


(58) 


(59) 


(60) 


(The  subscripts  are  really  necessary  only  on  the  interface.)   For 
any  function  f(x,y,t),   we  denote  the  values  on  the  interface  by 


f| 


def 


o 


f(x,  g(x,t),t) 


(61 


To  the  first  order  of  small  quantities,  the  equations  of 
fluid  dynamics  reduce  to 


^^±%^^      0  =  c2v20, 


(62) 


p  =  p 


o 


Po(^±"oll^' 


(63) 


where  the  upper  signs  are  to  be  used  for  y  **  g(x,t).  and  the 
lower  signs  for  y  <   g(x,t) 

The  boundary  conditions  at  the  interface,  to  the  first  order 
of  small  quantities,  are  first,   A«v, I   =  A«V  I   =  A'Vl  ,   where 
V  =  interface  velocity,  or 
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dt     o  dx  "^   '      dy  '  o  ' 

and  second,   ll  — §•  =  p,  I   -  P  I  ^   where  |j,  =  coefficient  of 
surface  tension,   or 

Ji_  ^^  (0     0  )  I  ^  p  (^  _  u   |-)0  I   -  p  (^  +  u  ^)0A 
2u    dxdy   '^H-    ^-  '  o   ^o  dt     o  dx  ^-'o    '  o  d  u     o  dx  ^+ '  o 
o     ^ 

■.65) 
[The  value  of  — §■  obtained  from  (64)  has  been  used  here.]   We  can 

ax^ 

eliminate   g(x,t)   from  (64)  to  give 


/O,      0\      !_/''-'  01    -t- 1 

^5t  +  %  ^'    W^o   -    ^^  ~   "^o     ^'   W~^o 


{66) 


Equations  (65)  and  {66)    serve  as  boundary  conditions  for  the  partial 
differential  equation  (62). 

The  solution  we  want  is  of  the  form 


ikx-t-cot 

(67 


0^   =  Re   A^  e^^-^^   cosh      ,<^(y   -   y^), 


ikx+cit 


0     =  Re   A      e^'^"^^''   cosh      k    (y  +-y^), 


o  ■ 


ikx-KJ3t 


g   =  Re      G   e"^-^\  (68) 


The  boundary  condition  that   a0/ay   should  vanish  at  the  rigid 
walls  at  y  =  +  y   is  satisfied  by  the  functions  in  (67).   The 
differential  equation  (62)  requires  that 
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(cD  +  Ik  u^)^  =  c^(«_^^-  k^), 


(o)  -  Ik  Uq)   =  c  (k_  -  k  ) 


(69) 


Upon  substituting  (67)  into  (65)  and  {66)    and  eliminating  A  , 
we  find  that 


coth  /c  y     ., 


coth  K,  y^    ., 

=  (CD  +  Ik  Uo)[-(a)  +  Ik  uo)  ^— ^  +  2^^]  • 

K+  qPo 


(70) 


Before  this  elimination,  equation  (65)  gives  the  ratio   A  /A   ,  viz., 


(00  +  ik  u^)K_A_sinh   «_  y^  =  -  (co  -  Ik  u^)  K_^A_|_sinh  «_^  y^;    ,   v 

lastly  G   is  obtained  by  solving  (64)  for  the  time-derivatives  and 
substituting  from  (67)  and  (68), 

-A.K.sinh  K ,  y^  +  A  k   sinh  k     y^ 
G  -  -i-i ±—^ ^^ =— ^  .      (72) 

2a) 

In  the  SUBMACH  code,  equations  (69)  and  (70)  are  solved  for 
CD,  K  ,   and  /c_   by  an  iterative  procedure  that  converges  rapidly 
if  u  /c   is  not  larger  than  about   3/4;   it  is: 
1)   initially,  take  ic^  =  K-   =    1^1, 
f^2)   solve  (70)  for  co  , 

5)   solve  (69)  for  .improved  values  of   |<^   and  (^_  , 
_4)   return  to  step  2). 
Subroutines  for  complex  arithmetic  are  used.   Equations  (69)  have 
two  solutions  for  i<  ,   equal  in  magnitude  and  opposite  in  sign. 
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and  two  similar  solutions  for  k    .   However,  only  even  functions 
of  K        and  of   (<:_   occur  in  the  whole  derivation,  so  the  choice 
of  sign  has  no  effect  on  the  solution.   There  are  also  two  choices 
in  the  solution  of  (70)  for  cd,   since  (70)   is  quadratic,  and 
the  two  corresponding  solutions  of  (69)  and  (70)  together,  after 
convergence,  can  he  obtained  from  each  other  by  replacing  cd  by 
-co  and  interchanging   k        and  ic_;    this  amor.nts  merely  to  a 
rotation  by  180   in  the  x-y  plane,  and  we  shall  therefore  be 
satisfied  with  one  of  the  solutions.   [Which  one  we  get  is  determined 
by  the  -  for  us,  unimportant  -  choice  made  by  the  subroutine  that 
takes  the  square-root  of  a  complex  number.] 

It  turns  out  that  cd   is  always  pure  real  or  pure  imaginary, 
but  no  use  was  made  of  this  in  the  coding. 

Although  the  primary  purpose  of  the  routines  described  in 
this  section  Is  to  provide  initial  and  comparison  data  for  the 
numierical  calculation,  they  also  give  direct  information  on  the 
stability  of  the  Interface.   Figures  10  and  11  shovj  how  od,      :;_    , 
and  i:_     vary  with  \i      for  a  typical  set  of  other  parameters. 

20,  Numerical  solution. 

The  calculation  utilizes  subroutines  of  the  MACH  code.   There 
are  two  regions  to  be  treated  by  regular  net  points,  a  slip  line, 
and  two  rigid  walls.   The  only  new  feature  is  the  periodicity  in 
x;   this  is  taken  care  of  by  having  extra  rows  of  net  points  (vertical 
rows,  in  Figure  9)  to  the  left  and  right  of  the  Interval   0  ^  x  ;^  1, 
and  supplying  them  with  values  of  u  v  p  p   from  corresponding  net 
points  inside  the  interval,  before  each  phase  of  the  calculation, 
after  which  the  standard  procedures  can  be  used  for  all  net  points 
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in  the  Interval,  extra  points  on  the  slip  line  are  similarly  provided. 
The  Instability  observed  In  the  MACH  code  was  also  present  here, 
as  expected.   The  code  was  then  modified  as  described  in  section  11, 
in  the  hope  that  that  method  would  be  stable.   Tests  so  far,  although 
very  brief,  seem  to  show  that  the  instability  is  still  present. 

PART  C   Proposed  Modifications  of  the  Method 

21 .  Introduction. 

I  now  plan  to  reformulate  the  problem  from  the  beginning, 
with  modifications  described  in  the  following  paragraphs,  and  program 
it  for  the  7090.   Parts  of  the  program  that  may  be  useful  in  a  wider 
class  of  fluid-dynamics  problems  will  be  put  in  a  subroutine  form 
convenient  for  future  tests  or  applications. 

22 .  Tentative  heuristic  analysis  of  the  instability. 

In  studies*  of  one -dimensional  shock  fitting  by  methods  analogous 
to  the  ones  used  here.  Dr.  John  Gary  has  encountered  a  phenomenon 
probably  related  to  our  instability.   In  one  dimension,  if   ^   denotes 
the  x-coordinate  of  the  shock  at  some  instant,  and  if  x.   denotes 
the  last  net  point  behind  the  shock,  the  derivative  of  a  function 
f(x)   is  estimated,  to  second-order  accuracy,  by  fitting  a  quadratic 
function  of  x   to  the  values  of   f  at   ^,  x.,   and  ^i_-i  -   ^^   ? 
and  X.   happen  to  be  very  close  together,  this  procedure  is  highly 
inaccurate,  since  a  small  error  of  the  value  of   f  at  either   | 
or  X.   can  make  the  graph  of   f(x)   bow  violently  upward  or  downward 
between  x,   and  x._-,    Gary  found  that  unless  measures  are  taken  to 
avoid  this  situation,  the  resulting  disturbance  seriously  impairs  the 

to  be  described  in  a  forthcoming  report  in  this  series. 


accuracy.   According  to  the  physical  interpretation  of  stability 
for  hyperbolic  difference  equations,  based  on  the  work  of  Courant, 
Priedrichs  and  Lewy  [7],  one  should  expect  trouble  whenever  the 
distance  between  two  points  used  in  the  computation  (such  as   ^ 
and  X.)   is  appreciably  less  than   cAt. 

J 

The  analogous  thing  in  two  dimensions  is  that  when  a  function 
f{x,j)      is  being  fitted  at  a  cluster  of  points  that  lie  close  to  a 
conic  section  in  the  x  -  y  plane,  small  errors  in  the  function 
values  can  cause  the  surface   z  =  f(x,y)   to  bow  violently  up  or 
down  in  the  region  inside  the  conic  section.   As  stated  in  section 
5,  this  is  not  such  a  very  rare  occurrence.   Presumably  trouble 
must  be  expected  whenever  there  is  a  conic  section  that  fits  the 
points  of  a  cluster  to  within  a  distance  appreciably  less  than   cAt. 

Neither  in  Gary's  work  nor  in  mine  did  this  disturbance  make 
the  shock  actually  unstable.   However,  a  shock  is  an  inherently 
stable  thing,  and  a  slip  line  is  not.   Hence,  I  surmise  that  the 
disturbance  may  be  the  cause  of  the  observed  instability.   In  any 
case,  it  is  clear  that  a  better  way  of  estimating  partial  derivatives 
is  needed. 

23 •  New  method  of  estimating  gradients. 

I  will  outline  a  procedure  that  is  more  direct  than  the  cluster- 
least-squares-quadratic-fit  method  and  in  which  too  closely  spaced 
points  are  avoided. 

As  a  first  step,  net  points  near  a  boundary  are  divided  into 
two  classes-   A  type   A  point  is  a  net  point  lying  closer  to  the 
boundary  than  a  certain  distance   5-,  ,   which  may  be  of  the  order 
•jj-  Ax   to  p-  Ax   --   the  difference  equations  are  not  used  at  such 
a  point;  a  type  B  point  is  a  net  point  at  least  one  of  whose  nearest 
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neighbors  lies  over  the  boundary  or  Is  of  type   A.   See  Figure  12. 

The  procedure  assumes  that  values  of  a  function  f  =  f(x,y) 
are  given  at  type  B  points,  at  ordinary  net  points,  and  at 
points  on  the  boundary  --  I.e.,  everywhere  except  at  type   A 
points c   The  partial  derivatives  of   f  will  be  obtained  at  type 
B  points  and  at  points  on  the  boundary <,   As  a  by-product,  interpo- 
lated values  of  f  will  be  obtained  at  the  type  A  points;  these 
will  be  needed,  later,  if,  owing  to  the  motion  of  the  boundary,  a 
type   A  point  becomes  a  type  B  point.   In  the  procedure,  no 
difference  quotient  is  ever  formed  for  a  pair  of  points  closer 
together  than  the  distance   5, ,   and  derivatives  are  always  estimated 
from  points  as  nearly  as  possible  in  a  straight  line;  the  formulas 
are  correct  to  C^((Ax)  ). 

The  first  step  is  to  supply  Interpolated  values  of  f  at  type 
A  points.   It  is  assumed  that  each  type   A  point  has  at  least  one 
neighbor  of  type  B.   [This  can  fail  to  be  so  in  a  narrow  region 
between  two  boundaries  that  meet  at  an  acute  angle  -  such  a  region 
must  receive  special  treatment.]   In  the  case  of  the  right-hand 
type   A  point  in  Figure  12,  the  type  B  point  below  it  is  used, 
and  the  final  Interpolation  is  with  r^espect  to  y.   If  the  vertical 
line  joining  these  two  points  is  extended  upward,  it  meets  the  point 
designated   X  on  the  boundary.   By  quadratic  interpolation  along 
the  boundary,  a  value  of   f   is  found  at  point  X.   If  this  same 
line  is  extended  downward,  it  usually  encounters  another  net  point 
in  the  same  region;  it  may,  however,  encounter  another  boundary 
first  (but  at  a  distance  ^  5.,  ),   where  a  value  of  f   can  be  obtained 
at  a  point  similar  to  X.   In  either  case,  we  now  have  three  points 
on  the  vertical  line,  and  we  can  obtain  f  at  the  type   A  point 
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by  quadratic  interpolation.   Other  cases  are  handled  similarly. 

In  the  second  step,  values  of  3—  and  3—  are  obtained 

ox       dy 

at  type  B  points.   Whenever  net  points  (including  type   A 
points)  on  opposite  sides  of  the  type  B  point  are  available,  an 
ordinary  centered  difference  quotient  Is  used.   Otherwise,  a  point 
similar  to  point  X  on  the  boundary  must  be  used  In  place  of  one 
(or  possibly  even  both)  of  the  neighbors,  and  then  the  approximation 
to  the  derivative,  must  take  the  unequal  spacing  of  the  points  Into 
account,  to  achieve  second  order  accuracy. 

In  the  third  step  (this  Is  actually  Independent  of  all  the 
other  steps  except  the  last),  values  of   A*»Vf  at  points  on  the 
boundary  are  obtained  by  differencing  with  respect  to  arc  length 
s   along  the  boundary;  if   A*,   the  unit  tangent  vector,  is  in 
the  direction  of  increasing   s,   we  have 


A*"Vf(x,y) 


=   ^f(x(s),y(s)) 


(73) 
Y 


where  Y   is  a  point  on  the  boundary.   The  finite-difference 
approximation  to   df/ds  must  take  possible  unequal  spacing  of 
points  along  the  boundary  into  account,  for  second  order  accuracy. 

In  the  fourth  step,  for  each  point  Y  of  the  boundary,  the 
nearest  point  of  type  B  is  chosen  (See  Figure  12);   if   s   denotes 
arc  length  along  the  straight  line  BY,   we  now  have  available  the 
values  of  f  at  B  and- at  Y,   and  the  value  of   df/ds  =  [x'Vf 
at  B,   where  |_l   is  a  unit  vector  along  BY.   This  suffices  to 
give  a  quadratic  fit  for  f(s)   and  from  this  the  value  of  ii.«Vf 
is  obtained  at  Y. 

In  the  last  step,  values  of   Sf/Sx,  and  bf/bj     at  Y  are 
obtained  as  linear  combinations  of   A*»Vf  and  |x«Vf   there. 
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If  the  angle  Z.A*,  \i      Is  too  small,  the  last  step  Is 
Inaccurate.   Hence,  we  make  an  exception  to  the  rule  of  the  fourth 
step  in  that  If  A*'1-l   Is  larger  than  a  certain  value,  when  the 
point  of  type  B  nearest  to  Y  is  chosen,  then  we  discard  this 
point  and  take  another  point  of  type  B. 

24.  Additional  Jump  condition  on  slip  line. 

In  section  15  it  was  pointed  out  that  although  the  rigid  wall 
might  have  been  treated  as  a  special  case  of  an  interface,  it  was 
instead  treated  by  a  reflection  principle.   Actually,  the  reflection 
principle  achieves  one  thing  that  interface  fitting  would  not  -  it 
makes  the  normal  component  of  the  pressure  gradient  A»Vp  vanish 
exactly  on  the  wall.   The  vanishing  of  A«Vp   is  an  exact  consequence, 
on  the  differential  equation  level,  of  equations  that  would  have  been 
used   if  interface  fitting  had  been  done,  but  it  is  presumably  only 
an  approximate  consequence,  on  the  difference  equation  level.   This 
may  be  of  importance,  since  the  calculation  appeared  to  be  stable 
near  the  walls. 

For  a  general  interface  (e.g.,  the  slip  line  in  the  Mach 
problem)  there  is  correspondingly  a  Jump  condition  for  A»Vp,   of 
which  no  use  is  made  in  the  fitting  procedure,  but  which  ought  to 
be  an  approximate  consequence  of  the  conditions  that  are  used.   It 
is  a  static  condition  (no  time  derivatives  appear)  and  intuition 
suggests  that  if  we  force  it  to  be  satisfied  exactly,  this  will 
strengthen  the  coupling  of  the  two  sides  of  the  interface  in  such 
a  way  as  to  improve  stability.   The  condition  is 
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[1  A.Vp]  -  [(A*.v)2]A*.;^  A  +  2[A^.v]  ;^(A.V)  ^  0,        (74) 


'Vf    o_«     'V'    U.O  'V'  .       ""u    -"u 


curvature         angular 
of  Inter-         velocity 
face  of  Inter- 

face 

where  the  symbol   [  ]   stands  for  the  jump  of  the  quantity  inside 

it.   In  the  first  term,  there  appears  the  analogue  of  the  quantity 

—  -v^  )      which  is  continuous  across  an  Interface  in  one-dimensional 
p  ox 

flow,  and  the  other  two  terms  represent  Jumps  of  centrifugal  and 
Coriolis  forces o 

Equation  (7^)  could  be  satisfied  by  a  slight  modification 
of  our  interface  fitting  procedure.   A  preliminary  step  of  that 
procedure  is  to  find  the  components  of   Vp   from  the  existing 
pressure  field  at  time  nAt   for  use  in  computing  various  quantities 
at  later  times.   Whether  the  numerical  differentiation  is  done  by 
the  method  of  section  5  or  by  that  of  section  23,  one  could  then 
add  a  small  amount  to  the  component  A»Vp   on  one  side  of  the  inter- 
face  and  subtract  the  same  amount  on  the  other  side,  to  satisfy  (7^). 

25.  Miscellaneous . 

In  the  next  formulation,  I  propose  to  use  the  Lax-Wendroff 
difference  equations  of  second  order,  and  to  apply  the  principle 
used  in  deriving  those  equations  to  all  the  fitting  procedures  as 
well.   This  may  not  be  strictly  necessary  in  the  Mach  reflection 
problem,  where  only  the  asymptotic  steady  state  is  desired,  but  it 
is  expected  that  the  methods  developed  here,  perhaps  even  the  actual 
subroutines,  will  be  used  also  in  other  problems. 

I  propose  to  treat  the  walls  by  Interface  fitting,  as  discussed 
in  sections  15  and  24,  in  order  to  reduce  the  number  of  special  methods 
that  have  to  be  programmed.   The  reflection  principle  would  not  be  of 
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much  use  for  curved  walls,  anyway. 

I  propose  to  allow  use  of  different  point  nets  for  different 
regions  of  flow.   In  the  Mach  reflection  problem,  the  region  under- 
neath the  slip  line  is  small  and  is  bounded  by  rather  sharp  corners, 
and  it  would  be  worthwhile  to  have  a  finer  net  there  than  in  the 
other  region. 

I  propose  to  use  the  method  of  section  23  for  estimating 
gradients.   A  way  of  treating  sharp  corners  still  has  to  be  worked 
out  in  this  method. 

I  propose  to  allow,  as  an  option,  the  enforcing  of  the 
additional  jump  condition,  described  in  section  24,  on  the  slip 
line  and  walls. 

Longer  range  objectives  should  include  better  treatment  of 
singular  points  in  the  x  -  y  plane.   There  are  six  such  points 
in  the  Mach  reflection  problem:  the  triple  point,  where  the  three 
shocks  meet,  2  points  where  a  shock  touches  a  rigid  wall,  a  point 
where  the  slip  line  touches  the  wall,  the  tip  of  the  wedge,  where 
some  kind  of  corner  flow  takes  place,  and  a  point  in  the  upper 
region,  so  far  ignored,  where  the  streamlines  converge  -  this  is 
a  singular  point  because  different  streamlines  have  different 
constant  values  of  entropy,  hence  of  pp~  ,   on  them. 

These  points  give  one -dimensional  singularities  in  the  three- 
dimensional  space-time.   There  is  surely  a  vast  variety  of  zero- 
dimensional  singularities;  in  fact   the  whole  Mach  reflection 
calculation  is  the  study  of  one  of  them,  for  the  event  that  occurred 
when  the  primary  shock  touched  the  tip  of  the  wedge  is  associated 
with  the  point   (0,0,0)   in  the  original   (X,Y,t)   system;  in  that 
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system,  It  Is  a  point  where  all  the  one-dimensional  singularities 
mentioned  above  come  together.   One  reason  for  trying  to  develop 
computing  methods  of  the  sort  here  described  is  to  study  other 
zero-dimensional  singularities » 

Another  longer  range  objective  is  to  give  better  Justification 
for  the  heuristic  arguments  about  the  stability  of  fitting  procedures 
for  shocks  and  interfaces » 
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TABLE  I  Summary  of  MACH  Calculation 


PHASE     SWEEPS 
SP,  SPA 

SPA       b'dles 

topology 

SP       net 

►  IC        net 

IC 
ID 
ID 

2 

SPA 

6A 

6B 

6C 

6d 

7        net,  b'dles 
output  A  net,  b'dles 

output  B  net,  b'dles 
streamlines 
J>A  net 


b'dles 

long 

list 

long 

list 

long 

list 

net. 

often 

net 

b'dles 

shock 

slip 

line 

QUANTITIES  CALCULATED,  STEPS 

general  constants,  3-shock  conflg.,  tables  of 
functions  used  In  shock  fitting. 

Initial  data  on  shock  and  slip  linCo 

region  indices. 

initial  data  on  net,  long  list  entries  for 
permanent  clusters. 

net  point  tags,  remaining  long  list  entries, 
^n+1/2  _^     ^n+3/2^ 

ident.  of  net  point  nearest  to  each  bdy  pt. 

cluster  descriptions, 

interpolated  new  data  for  any  net  point  that 
has  changed  region. 

coefficients  of  least-squares  fit. 

(n  =  0  only)  smooth  values  of  u,  v  by  relaxation. 


partial  derivatives  of  u,  v,  p, 

n+1 
p     on  neto 

T.n+1 


and 


P 


n+1 


P    -  surf.  tens. 


.n+1 


n+1 


(A-v) 


n+1 


c     on  shock. 


surf,  tens,  p''^^,  p''^^  (A.v)'^+^,  c 


^  P 

n+1 


n+1 


on  slip  line, 


n+1 


n+1 


3B 


b'dles 


p    -^     p,    p    —^  p   on  net;   n  +  1  ->  n. 
short  print-out,  picture  (skip  output  B  if 
switch  6  is  up) o 

longer  print-out,  entropy  test,  Bernoulli  test, 
picture  of  streamlines. 

partial  derivatives  of  p"^  and  p"^,  u^^  /^, 

'^.n+J  /2 

V       on  net, 

Jn+1/2^ 
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PHASE 

SWEEPS 

3C 

shock 

3D 

slip  line 

4A 

net 

4b 

b 'dies 

4c 

shock 

4d 

slip  line 

5A 

net 

5B 

b'dies 

5C 

shock 

5D 

slip  line 

lA 

b 'dies 

IB 

b 'dies 

topology 

QUANTITIES  CALCULATED,  STEPS 

'^.n+1/2  ~n+l/2      ,   , 
u   '     ,    V        '         on  shock. 

'N^n+1/2  ~n+l/2 

u   '     ,    V        '         on  slip  line" 

partial  derivatives  of  u  and  v,  p   ,  p 

~n+l  -vn+l      ,   , 
p    ,  p     on  shock. 

p    ,  p     on  slip  line. 

n+1/2        n+1/2  , 

u  ,    V        '         on  net. 

pn-fl/2^ 

n+1/2    n+1/2   ,  „   n+1/2      ,   , 
u      ,  V   '     ,    A»v,  p    '         on  shock. 

same  on  slip  line. 

n+1  /2 
new  positions  of  shock  and  slip  line;  v   '      — >  v 


a_» 


new  normal  vectors   A. 

region  indices,  table  of  u,  v,  p,  p,  x,  y   if 
switch  6  down  (printed) . 

IE        b'dies         interpolated  values  at  new  b ' dy  points,  if  any 

are  needed. 

—  return  to  phase   IC  — 
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Fig.  1 
Mach  reflection  In  a  shock  tube. 
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Pig.  2 
Net  points  and  boundary  points 
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Fig.  3 
Point  clusters  used  In  the  least-squares  quadratic  fitting 
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Fig.  4 
A  characteristic  plane  tangent  to  the  sonic  cone 
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Fig.  5 
The  initial  configuration, 


slip  line 


clusters 


wedge 


Fig.  6 
Typical  clusters  at  the  surface  of  the  wedge 
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Fig.  7 
A  special  cluster  at  the  surface  of  the  wedge 
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Fig.  9 
Diagram  of  the  SUBMACH  problem. 
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Fig.  12 
Type  A  and  Type  B  net  points . 
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